Compressive Optical Deflectometric Tomography: 
A Constrained Total- Variation Minimization Approach. 
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Abstract 

Optical Deflectometric Tomography (ODT) provides an accurate characterization of 
transparent materials whose complex surfaces present a real challenge for manufacture and 
control. In ODT, the refractive index map (RIM) of a transparent object is reconstructed by 
measuring light deflection under multiple orientations. We show that this imaging modality 
t-H can be made compressive, i.e., a correct RIM reconstruction is achievable with far less ob- 

servations than required by traditional Filtered Back Projection (FBP) methods. Assuming 
a cartoon-shape RIM model, this reconstruction is driven by minimizing the map Total- 
Variation under a fidelity constraint with the available observations. Moreover, two other 
realistic assumptions are added to improve the stability of our approach: the map positiv- 
ity and a frontier condition. Numerically, our method relies on an accurate ODT sensing 
model and on a primal-dual minimization scheme, including easily the sensing operator and 
{vq the proposed RIM constraints. We conclude this paper by demonstrating the power of our 

t-H method on synthetic and experimental data under various compressive scenarios. In par- 

ticular, the compressiveness of the stabilized ODT problem is demonstrated by observing a 
typical gain of 20 dB compared to FBP at only 5% of 360 incident light angles for moderately 
noisy sensing. 



Keywords: Optical Deflectometric Tomography, Refractive index map, Total Variation Min- 
imization, NFFT, Primal-Dual Optimization. 



1 INTRODUCTION 

Optical Deflectometric Tomography (ODT) is an imaging modality that aims at reconstructing 
the spatial distribution of the refractive index of a transparent object from the deviation of 
the light passing through the object. By reconstructing the refractive index map (RIM) we are 
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able to optically characterize transparent materials like optical fibers or multifocal intra-ocular 
lenses (IOL), which is of great importance in manufacture and control processes. 

ODT is attractive for its high sensitivity and effective detection of local details and object 
contours. Also, ODT is insensitive to vibrations and it does not require coherent light. More- 
over, compared to interferometry, ODT can be applied to objects with higher refractive index 
difference with respect to the surrounding solution. 

The deflectometer used for the data acquisition is based on the phase-shifting Schlieren 
method [17] . For each orientation of the sample, two-dimensional (2-D) mappings of local light 
deviations are measured. As these light deviations are proportional to the transverse gradient 
of the RIM integrated along the light ray, ODT is able to reconstruct the RIM from the angle 
measurements. 

First works on deflectometric tomography [3 16 have focused on the use of common re- 
construction techniques like the Filtered Back Projection (FBP). They have proved that FBP 
provides an accurate estimation of the RIM when sufficient amount of object orientations is 
available. However, when we consider a scenario with a limited number of light incident an- 
gles, and in the presence of noise in ODT observations, FBP induces the apparition of spurious 
artifacts in the estimated RIM, lowering the reconstruction quality. 

Inspired by the Compressed Sensing (CS) paradigm (6 12 , which demonstrates that few 
measurements are enough for an accurate reconstruction of sparse signals, recent works in 
ODT have started to exploit sparsity to reconstruct the RIM from few acquisitions, i.e., in the 
presence of an ill-posed inverse problem. Foumouo et al. [XT] and Antoine et al. [I] have used a 
sparse representation in a B-splines basis and regularized the problem using the ^i-norm. The 
reconstruction was performed using iterative schemes and the results show that, although the 
method is capable of reproducing the shape of spatially localized objets, the image dynamics is 
not well preserved and the borders are smoothened. 

Although sparsity based methods are new in ODT, they have been used in other applications, 
such as Magnetic Resonance Imaging (MRI) [25], Absorption Tomography (AT) |30,32|, Radio 
Interferometry [35] and Phase-Contrast Tomography 10 IT], for the image reconstruction when 
few amount of acquisitions is available. More details are given in Sec. [4] 

An additional problem that rises in all physical applications is the estimation of the sensing 
operator that fits better the physical acquisition. Most operators present a non ideal behavior, 
which conditions the numerical methods to solve the inverse problem. In tomographic problems, 
this operator requires to map spatial data in a Cartesian grid to Fourier data in a Polar grid. 
Previous works have used gridding techniques to interpolate data from a polar to a cartesian or 



pseudo polar grid before applying the Fourier Transform 19 36 . However, the error introduced 



when using these techniques is not bounded and introduces an uncontrolled distortion. 



1.1 Contribution 

In this work, we show that ODT can be compressive and robust to Gaussian noise. We propose a 
constrained method based on the minimization of the Total- Variation (TV) norm that provides 
high quality reconstruction results even when few acquisitions are available and in the presence 
of high level of Gaussian noise. 

We assume the RIM is composed by slowly varying areas separated by sharp transitions 
corresponding to material interfaces. Such a behavior is known to be represented by spatial 
functions having a small TV norm. Other sparsity basis, such as the B-splines basis in [T]|17|, 
tend to be complicated to use and do not provide a complete characterization of the RIM, while 
the TV norm represents a more simple and accurate model of the actual RIM. 

Therefore, we propose the use of the TV norm to regularize the ill-posed inverse problem of 
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estimating the RIM from low amount of noisy data. To account for the noise and the raw data 
consistency, we add an £2 data fidelity term adapted to Gaussian and uniformly distributed 
noise. 

The proposed method offers the flexibility to work with more than one constraint. In this 
work, we add two more constraints based on some general prior information on the RIM in 
order to converge to an optimal solution: (i) the RIM is positive everywhere and (ii) the object 
is completely contained in the field of view of the ODT experiment. These added constraints 
are proved to provide an unique solution to our ODT problem. 

As for the operator, we use the NFFT algorithm^ a fast computation of the non-equispaced 
DFT. This algorithm provides an efficient estimation of the polar Fourier transform with a 
controlled distortion regarding the true polar transform. 

The compressive ODT problem is solved by means of the primal-dual algorithm proposed 
by Chambolle et al. (8j. The results are compared with a common FBP approach, showing the 
proposed method outperforms FBP in terms of compressiveness, noise robustness and recon- 
struction quality. 



1.2 Outline 

In Sec. [2] we provide a brief background on optical deflectometric tomography, describing also 
the experimental setup used for the data acquisition. Then, the ODT discrete model is pre- 
sented in Sec.|3j In Sec. [4] we depict related works on tomographic reconstruction, which provide 
a basis on the methods adopted to recover the RIM: the commonly used FBP method and a 
proposed regularized method that we call TV-^2- These methods are described in Sec. [5| In 
Sec. [6] we present the identification and estimation of the noise in both synthetic and experimen- 
tal data, and the analysis of the noise impact when comparing common absorption tomography 
and deflection tomography. Sec. [7] presents the numerical methods we use to recover the RIM 
from the noisy measurements by means of the proposed regularized formulation. Finally, in 
Sec. [8] some reconstruction results are shown, focusing first on the comparison between common 
tomographic and ODT reconstructions, and then on the comparison of the reconstruction meth- 
ods on the basis of quality reconstruction and convergence for both synthetic and experimental 
data. 



1.3 Conventions 

Most of domain dimensions, e.g., M, N, N\, N%, . . ., are denoted by capital roman letters. 
Vectors and matrices are associated to bold symbols while lowercase light letters are asso- 
ciated to scalar values. To give an example, for some dimensions D,D' G No, we write 
u = (u\, ■ ■ ■ ,ud) t G M. D (using the transposition or we define the matrix <& G M. D xD 

with matrix entry <!>,.,■ G R. Therefore, the i th component of a vector u is denoted either u\ or 
(u)i, the later notation being useful if u brings some extra notations (e.g., us), while U{ will 
denote for instance the i th vector of a set of vectors 112, ■ ■ ■ }. The identity matrix in R D is 
denoted 1 D , or simply I when its dimension is clear from the context. The set of indices in M. D 
is [D] = {1, • • • , D}. Scalar product between two vectors u, v G M. D for some dimension D G N 
is denoted equivalently by u T v = u ■ v = (u,v). For any p ^ 1, || ■ || p represents the £ p -novm 
such that \\u\\p = ~Yli\ui\ p with || • || = || • ||2- The notation \\u\\q = #suppu represents the 
cardinality of the support supp u = {i : Uj 7^ 0} C [D]. For a subset S C [D], given u G ~SP and 
$ G MP xD , us G (or $>s) denotes the vector (resp. the matrix) obtained by retaining the 
components (resp. columns) of u (resp. 3>) belonging to S. Alternatively, we have us = Rsu 



This toolbox is freely available here http://www-user.tu-chemnitz.de/~potts/nfft/ 
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Figure 1: (top) The Deflectometric model, (bottom) A scheme of a Phase-Shifting Schlieren Defiectometer (in 
refraction mode) measuring light ray deflection angles by encoding light deviation into light intensity variations. 



or = where R$ := (I,s) T G {0, l}# SxD is the restriction operator. The kernel of a 

matrix $ is defined as ker(<&) := {x G MP : &x = 0, x ^ 0}. The norm of an operator K is 
defined as \\K\\ = max{||i^a;|| : x G with ||a;|| = 1}. We denote r°(V) the class of proper, 
convex and lower-semicontinuous functions from a finite dimensional vector space V (e.g., M. D ) 
to (— oo,+oo] [9 15 . We denote by (x) + the non-negativity thresholding function, which is 



defined componentwise as (xi)+ = (xj + |xj|)/2. We denote by ic{x) the indicator function of 
the set C, which is equal to if x G C and +00 otherwise. We denote by int Q the interior of a 
set f2, which consists of all points of £1 that do not belong to the boundary 50 of O. 



2 OPTICAL DEFLECTOMETRIC TOMOGRAPHY 

In this section, the principles of optical deflectometric tomography are explained and completed 
with a brief description of the experimental setup used for actual deflectometric acquisition. 



2.1 Principles 

Optical Deflectometric Tomography aims at infering the refractive index spatial distribution 
(or refractive index map - RIM) of a transparent object. This is achieved by measuring, under 
various incident angles, the deflection angles of multiple parallel light rays when passing through 
this transparent object (see Fig. [TJ-(top)). The (indirect) observation of the RIM, allowing its 
further reconstruction, is guaranteed by the relation between the total light ray deflection angle 
and the integration of the transverse gradient of the RIM along the light ray path |3] . 

In this work, we restrict ourselves to two-dimensional ODT by assuming that the refractive 
index n of the observed object is constant along the e3-axis for a convenient coordinate system 
{ei,e2,es}, i.e., dn(r)/dr% = (with r = (ri,r2,r^) T G M 3 ) and deflections occur in the ei-e2 
plane. This strong restriction can be softened by assuming that the variations of n along the 
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e% axis are slow compared to those occurring in the {ei, domain. 

Consequently, the direction e% can be forgotten. Given the refractive index n : r = 
(ri,r2) T G M 2 — > n(r), for a particular light ray trajectory 1Z = {r(s) : s G M} C M 2 
parametrized by r(s) = (ri(s), r2(s)) T G M 2 with s describing its curvilinear parametei[^J the 



deflection angle is provided by [22] 

A(K;n) = ± [ [Vn{r(s)) • p( s )] ds, (1) 

where n r is the (constant) reference refractive index of the surrounding medium, p(s) is the unit 
normal vector to the curve 1Z in r(s) such that p(s) ■ ^r(s) = 0, and the quantity Vn(r(s)) -p(s) 
is actually the transverse gradient of n along 1Z (the directional derivative of n along p) . 

Of course the light ray trajectory 1Z depends itself on n as computed by the light ray equation 
2^(n^r(s)) = Vn established from the eikonal equation [Zj. However, for small deflection 
angles A, we can adopt the (first order) paraxial approximation and assume the trajectory 1Z 
is a straight line. 

In this simplified model, the 2-D RIM is measured by the deflection angle A(r, 9;n) of a 
light ray TZ(8,t) = {r G M? : r • p 9 = r}, where r £ 1 is the affine parameter determining 
the distance between 1Z and the origin, 6 G [0, 2tt) is the incident angle of 1Z with the e\ 
axis, and p e = (— sin 9, cos 9) T is the (constant) perpendicular vector to the light ray direction 
tg = (cos9,s'm9) T . 

Consequently, removing the explicit trajectory dependence in A, the simplified ODT model 
depicted in Fig. [IJ (top) reads 

A(r,9):=l [ Vn(r Tie (s))-p e ds 
Jr 

= i I (Vn(r)-p,) 5(T-r.p e )d 2 r, (2) 

where r Tt g(s) = stg + rp g G 1Z, and the Dirac distribution turns the line integral into an 
integration over M 2 . In short, the above equation relates the deflection angle A(r, 9) to the 
Radon transform of the transverse gradient of n within the paraxial approximation. 

As for traditional parallel absorption tomography 23 1 , there exists a Deflectometric Fourier 
Slice Theorem (DFST) that relates the 1-D (radial) Fourier transform of the deflection angle 
along the affine parameter r, i.e., 

y(u,0) := / A(r,#)e- 2 -™dT, (3) 

to the 2-D Fourier transform of the RIM. Mathematically, the DFST establishes the following 
equivalence [19] (proved in Appendix |A|) : 

y(u,,6) = 2 -^n{u Pe ), (4) 

where fi(fc) = J" R2 n ( r ) e~ 2mk ' r d 2 r stands for the 2-D Fourier transform of n. Hereafter, the 
value ui is called affine frequency. 

We may remark from Q that there are different ways to cover the 2-D frequency plane with 
deflectometric measurements. This can be observed from the following symmetry relations for 

3 By a slight abuse of notation, r denotes any points in R 2 while r(s) represents a particular curve in R 2 
parametrized by s £ R. 
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n G M, r,w e M, 61 e [0,2tt): 

A(t, (fl + vr) mod 2vr) = -A(-t,6>), (5a) 

y(w,(6> + vr) mod 2vr) = 9), (5b) 

y(u,0)=y*(-u,,e), (5c) 

where the two first relations come from the change Pq +7T = —pg in ^ and Q, and the last 
equation is the well known Fourier conjugate symmetry for real spatial functions. 



In particular, from the symmetry (5c), we can restrict uj to positive values by taking in the 



whole circle [0, 2tt). Or alternatively, from (5b) and (5c), we can deduce y(uj, 9) = —y*(uj, (9 + ir) 
mod 2tt) and restrict 9 to the half circle [0, tt) with uj G R. We insist on the fact that this 
symmetry is only preserved when the paraxial approximation is validated. 

When comparing the relation Q with the standard tomographic Fourier Slice theorem 
(FST) 123], the main difference is provided by the transverse gradient in the deflectometric 
relation Q, which results in multiplying by 27riu>/n r the RIM Fourier transform. In particular, 
from ([2]) or from Q (since uj vanishes on the frequency origin) we see that the ODT sensing is 
blind to constant RIM. As we will see in Sec. [5j this has an impact on the formulation of the 
RIM reconstruction since the RIM can only be estimated relatively to the reference RIM n r . 

Let us conclude this section by insisting on the impact of the equivalence Q. Similarly 
to the use of the Fourier Slice theorem in common (absorption) tomography, Q is of great 
importance for defining a discrete ODT sensing model which can be computed efficiently in the 
Fourier domain given a discretized refractive index map n. 

2.2 Deflection Measurements 

Experimentally, the deflection angles A can be measured by phase-shifting Schlieren deflection 
tomography (schematically represented in Fig. [I]- (bottom)). We briefly explain this system for 
the sake of completeness in order to set the experimental background surrounding the actual 
deflection measurement process. More details can be found in (3}[T7|[l9j[2l] . 

This system proceeds by encoding light ray deflection a in intensity variations. A transparent 
object is illuminated with an incoherent uniform light source Io modulated by a sinusoidal 
pattern m using a Spatial Light Modulator (SLM). From classical optics, the light deviation 
angle a is related to a phase shift Ax = / tana, where / is the focal length of Lens 1. This 
phase shift is associated with the intensity variation thanks to the modulation m as I(—t, 9) = 
m(Ax)Io. These intensity variations are processed by phase shifting methods for recovering 
the deflection measurements a = A(t,9) for each couple of parameters (r, 9) E R x [0, 2n]. 
Up to some linear coordinate rescaling, the affine parameter r corresponds to the horizontal 
pixel coordinate in the 2-D CCD detector collecting light (assuming the object refractive index 
constant along the CCD vertical direction). This correspondence is implicitly allowed by the 
telecentric system formed by the combination of Lens 2 and 3. The pinhole guarantees that 
only parallel light rays outgoing from the object are collected. Rather than rotating the whole 
incident light beam around the object, it is this one which is rotated by an angle —9 along 
an axis parallel to the CCD pixel vertical direction [3J. Finally, since the system is invariant 
under time inversion, i.e., under light progression inversion, measuring the deflection angle a 
in Fig. [IJ- (bottom) is equivalent to measuring the same angle in Fig. [TJ-(top). 

3 DISCRETE FORWARD MODEL 

In order to reconstruct efficiently the RIM from ODT measurements, recorded data must be 
treated appropriately considering, jointly, the data discretization, the polar geometry of the 
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ODT sensing and unavoidable measurement noises. In this section, the discrete formulation 
of the ODT sensing and the construction of the forward model from the recorded data are 
presented. 



3.1 Discrete domains 

Let us first assume that the object of interest is fully contained in a square field-of-view (or 
FoV) Cl C M 2 centered on the spatial origin. The physical dimensions of this FoV can be 
provided by the Deflectometric device itself. In other words, the RIM is constant and equal to 
the reference index n r outside of 0. This involves also that the deflection measurement vanishes, 
i.e., A(t, cj) = 0, if |r| is bigger than the typical width of ft, in a section of direction 9 + tt/2. 

We can consider a spatial sampling of as follows. We define a Ao x Nq 2-D Cartesian grid 
of N := Nq pixels as 

Cn = {r m , n ■= (mSr,n5r) : -N /2 ^m,n< N /2}, 

where the spatial spacing 5r is adjusted to ft and to the resolution by imposing Q = [— ^Nq Sr, |iVo Sr] x 
[-±N 5r,lN 6r]. 

Second, as the deflectometric experiments provide evenly sampled variables r and 9, A is 
measured on a (signed) regular polar coordinate grid 

V M := {{r s ,e t ) : -(N T /2) ^ s < {N T /2), 0^t< N }, t s := s8t, 9 t :=tS9, 

of size M := N T Ng, with N T the number of parallel light rays passing through the object 
(assumed even), Ng the number of incident angles in ODT sensing, 5t and 59 = tt/Nq the 
distance between two consecutive affine parameters and angles respectiveljj^J 

The value 6t can be known experimentally from the pixel size of the CCD detector in a 



Schlieren Deflectometer (Sec. 2.2). Moreover, the value 5t and the resolution iVV are also related 
to the FoV ft so that 5tN t ~ 5rNo. Since there is no reason to ask more resolution to the 
sampling Cn than the available in the affine variations of the ODT measurements, we will work 
with 5t ~ 5r and N T ~ iVo- 

Third, in this discretized setting, the affine frequency w in Q must also be sampled with 
iV T values. As described in the next section, this comes from the replacement of the (radial) 
Fourier transform in ^ by its discrete counterpart. This leads to a (signed) frequency polar 
grid of same size 

V M := {(ov,0t) = ~(N T /2) < s' < (N T /2), 0^t< Ng}, u s > := s' 8oo, 9 t :=tS9, 

with 5uj = 1/(N t 5t). As it will become clearer in the following, only half of this polar grid will 
be necessary to bring independent observations of the RIM, i.e., we will often work on 

:= {(uv,0 t ) : s: s' < (N T /2), 0^t< Ng}, 

with #P+ = M/2. 



3.2 Discretized Functions 

From the discrete domains defined above, the continuous RIM n observed in the experimental 
FoV is discretized into a set of N values n(r mjn ) from the coordinates r m) „ G Cn- This descrip- 
tion can always be arranged into a one dimensional N- length vector n € WL N , given a convenient 

4 Notice that 80 is not set to 2n/Ne since r is allowed to be negative. 
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mapping between the components indices of n = (ni, • • • , njv) T and the pixel coordinates in C/y. 
This "vectorization" provides us a simplified representation of any linear operator acting on the 
sampled RIM n as a matrix multiplying n. 

For the different functions discretized on Vm or on Pjf, we use the same vectorization trick, 
namely, a function u defined on Cn and sampled on Vm is associated to a vector u G R M with 
the right correspondence between the components of u and the polar coordinates in Vm (and 
similarly for function defined in the Fourier domain and sampled on Vm)- 

Therefore, the ODT observations { A(r, 9) : (r, 9) G Vm} are gathered in a vector z G R M 
with Zj = A(r s ,9t) for a certain index mapping j = j(s, t) G [M]. In this discrete representation, 
the equivalent transformation of ^ reads 

y comp = (^r5T)F™ d z, (6) 

where y corap G C M is associated to a (vectorized) sampling of y on 

V M , and F md : R M -> C A/ 

performs a 1-D DFT on the radial r-variations of z, i.e., 



(F rad z) fc(sV) = -±= z j(s,t) e 



2wiss'/N T 



for a vectorized index k = k(s',t). In other words, if 5t is sufficiently small, e.g., if A(t, 0) is 
band-limited with a cut-off frequency smaller than 1/Sr = N t 5uj, then 

y«s>,t) = ^§r E A (^ 1 )^ l * ),s7< * ) 

= ^A(r s ,^)e- 2 ^' 5r 

s 

~ y(u s >,9 t ), 

using a Riemann sum approximation of ([3| and knowing that A vanishes outside of 0. 

/comp 



Despite the fact that y comp belongs C M ~ M 2M , this vector brings only M independent real 



observations of n G M. N . This is due to the central symmetry ( 5c ) induced by the realness of z 
and which allows us to consider y corap only on V M . 

This is clarified by the definition of the useful operator : C M — > M M which perform the 
two following linear operations. First, it restricts any vector £ G C M to the indices associated 
to the half grid V M . Second, it appends the M/2 imaginary values of the restricted vector 
to its M/2 real values in order to form a M -length vector in M M . The adjoint operation 0*, 
which is also the inverse of for vectors in C M respecting the Hermitian symmetry, is obtained 
easily by first reforming a M/2- length complex vector from the separated real and imaginary 

?+ 

M 



parts, and by inserting the results in C M according to the indices of V M and by completing the 



information in Vm \ Pm with the central symmetry (5c) 



Consequently, thanks to 0, we can form the real vector 

y = &y comp = (y r N~r5r)(@F^ d )z G R M , (7) 

with (0F rad ) : R M -> R M . We call y the Frequency Deflectometric Measurements (FDM) 
vector and we will most often use it as our direct source of ODT observations instead of z. 

3.3 Forward Model 

We can now explain how we use the DFST relation Q for defining the forward model that links 
any discrete 2-D RIM representation to its FDM vector. 



S 



In a previous work 19 , the data available in the frequency polar grid Vm were first inter- 
polated to a Cartesian frequency grid in order to reconstruct the 2-D RIM using a Discrete 
Fourier Transform (through the FFT algorithm). However, the polar to Cartesian frequency 
interpolation introduced a hardly controlled colored distortion. 

We prefer here the use of another operator F : ~R N — > C M performing a non-equispaced 
Discrete Fourier Transform (NDFT) for directly relating functions sampled on the Cartesian 
spatial grid Cn to those sampled on the polar frequency grid Vm- 

More precisely, given a function / : Vt — > C sampled on Cat, the NDFT0 computes 

JVo-lJVb-1 

/>)= E E f(r m ,n)e- 2mk - r ^, (8) 

m=0 n=0 

on the M nodes k of Vm ■ Gathering all the values of / and / into vectors / G C M and / G M. M 
respectively, this relation is conveniently summarized in the following equation 

f = Ff, with Fij = e -** ihi - r i. (9) 

where the matrix F G C Mx stands for the linear NDFT operation. Its internal entry in- 
dexing follows the one of the components of / and /. We explain in Appendix [b] how the 
Non-equispaced Fast Fourier Transform (NFFT) algorithm allows us to compute efficiently in 
0(Nlog(N/e)) the multiplications Fu and F*v for any u G l w and v G C M , with a controlled 
distortion e with respect to the true NDFT. 

The action of F on a discretized RIM n is related to the continuous Fourier transform of n 
as follows. Let j = j(s',t) G [M] be the j th point kj of the (vectorized) grid Vm associated to 
the polar coordinates (ui s r,9t). Then, for a sufficiently small Sr, 

N -1N -1 

( F ")Ks>,t) = E E <r m<n )e~ 2mk r^ « ^ n(avp flt ). (10) 

m=0 n=0 

To take into account the multiplication by 2irioj/n v in Q and the existence of a factor 



l/(8r) in the equivalence (10), we introduce the diagonal operator D G R MxM defined as 



D = 27T 'f r r) diag(w (1 ), • • • ,W(Af)), 

where corj\ refers to the w-coordinate of the j th point of Vm, if as before j = j(s',t) is 
associated to (uv, Ot) G Vm then uju\ = uj s >. The operator D models the effect of the transverse 
gradient in the Fourier domain. 

In parallel to the discussion ending Sec. |3.2[ we a lso restrict the action of DF to the domain 



V M . Consequently, using the operator (Sec. 3.2), the final linear forward model linking the 



real FDM to the 2-D NDFT of the discrete RIM n reads 

y = (@DF)n + r] el A/ . (11) 

The additional noise rj G M M integrates the different distortions induced by the numerical 
computations (e.g., the NFFT inaccuracy when estimating the NDFT F, see Appendix [B| ), 
the model discretization (e.g., the different Riemann approximations), the discrepancy to the 
paraxial approximation Q, and the actual noise corrupting the observation of z. A detailed 
noise estimation is provided in Sec. [6} 



This NDFT formulation is strictly equivalent to the one given in 24 where St = Sr = 1 
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Notice that in the absence of noise (r) 



0), the model (11) could be easily turned into 

Indeed, forgetting 



a classical tomographic model where the DC frequency is not observed, 
the transparent action of 0, if the frequency origin has been carefully removed from Pm-, 



then D is invertible with D 
tomographic problem 



y 



diag(w (1 J, 



(M) 



D 



y 



Fn. 



), and we can solve the common 

(12) 



However, as we present in Sec. 8.1.1 this transformation is not suited to noisy ODT sensing. 
Even for a simple additive white Gaussian noise r) (or AWGN), the multiplication by D~ l 
breaks the homoscedasticity of rj, i.e., the variance of each (D~ 1 r])j varies with j £ [M]. This 
interferes with common reconstruction techniques used in classical tomography. Obviously, a 
noise whitening could be realized for stabilizing such methods but at least for AWGN, this 



strictly amounts to solve directly the model (11) 



4 RELATED WORKS 

This section describes some recently proposed methods for tomographic reconstruction in the 
domains of differential phase-contrast tomography and common absorption tomography. 

In the domain of differential phase-contrast tomography, we face the problem of recovering 
the refractive index distribution from phase-shifts measurements. These measurements are 
composed by the derivative of the refractive index map, inducing the apparition of the affine 
frequency oj when using the FST, as it happens in the ODT sensing model described in Sec. [2] 
In this domain, Pfeiffer et al. |28| have used the FBP algorithm to reconstruct the refractive 
index map from a fully covered set of projections. Cong et al. [10[|11| have used different 
iterative schemes based on the minimization of the TV norm to reconstruct the refractive index 
distribution over a region of interest. These methods are accurate and provide similar results, 
but the iterative scheme based on the TV norm has proved to be better than FBP when the 
amount of acquisitions decreases. 

In common Absorption Tomography (AT) we deal with the reconstruction of the absorption 
index distribution from intensity measurements. As these measurements are directly related to 
the absorption index, the AT sensing model does not include the affine frequency oj. In this 
domain, several works have exploited sparsity based methods. Most recent works in AT have 
focused on promoting a small TV norm [301321. Sidky et al. [32] use a Lagrangian formulation for 
the tomographic reconstruction problem, promoting a small TV norm under a Kullback-Leiber 
data divergence and a positivity constraint. They aim at reconstructing a breast phantom from 
60 projections with Poisson distributed noise. For this, they use the primal-dual optimization 
algorithm proposed by Chambolle et al. [8]. The method results in high quality reconstruction 
compared to FBP but with a convergence result that is highly dependent on the Lagrangian 
parameter choosen. 

Ritschl et al. 1 30 1 use a constrained optimization formulation to reconstruct the absorption 
index from low amount of clinical data in the presence of metal implants and Gaussian noise. 
This problem is solved by means of an alternating method that allows then optimizing separately 
the raw data consistency function and the sparsity cost function, without the need of prior 
information on the observations. The fast convergence of the method is based on the estimation 
of the optimization steps. The gradient descent method is used to minimize the TV norm and 
the consistency term is minimized via an algebraic reconstruction technique. The method is 
proven to give better results than FBP. 

These works have proved that some tomographic applications can be made compressive in 
the CS sense (6, 12 . They show that it is possible to perform an accurate tomographic image 
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reconstruction from a number of samples that is smaller than the desired image resolution, when 
the image is sparse in some basis, i.e., the image expansion in that basis contains only a small 
number of nonzero coefficients. However, it is important to notice that most of these works have 
considered Cartesian grids, an approximation of the actual sensing in a polar grid, and plus 
they attack different problems than ODT, where the sensing model and type of measurement 
change from one to the other. One important difference between AT and ODT sensing models, 
relies on the presence of the affine frequency as materialized by a diagonal operator D; whose 



impact is analyzed in detail in Sec. 6.3 and Sec. 8.1.1 for perfect and noisy sensing. 



5 REFRACTIVE INDEX MAP RECONSTRUCTION 

Two methods can be considered for recovering the discrete RIM n 6 R N . The first is the Least 
Square solution, also known as Filtered Back Projection (FBP), and the second is a regularized 
approach called TV-^2 minimization. Both procedures are described in detail hereafter. 

Least Square Reconstruction 

Given any linear sensing model 

y = $x € R M (13) 

of some image (vector) x E M. N by a sensing operator 3? £ M. MxN with M ^ N, the Least 
Square solution x\ s is obtained by solving 

a; is = argmin ||u||2 s.t. y = &u. (14) 



This minimization amounts to regularizing the general inverse problem aiming at reconstructing 



x from ( 13 ) by promoting a minimal ^2-norm of the solution. As a matter of fact, X\ s is obtained 
in a closed form by applying the Moore-Penrose pseudoinverse of $ on the observed data y, 
i.e., 

x\ s = &y = $*($$*) _1 y. 

As explained in the previous section, the ODT sensing operator is $ = &DF, the ODT 
sensing model is y = &DFn and we denote by ni s the corresponding LS solution. Notice 
that if our forward model followed a classical AT sensing rather than an ODT one, it would be 
reduced to a simple partial Fourier sampling of n by $ = &F, without the action of D. 

Notice that the LS solution is equivalent to a discretization of the well known (analytic) 
Filtered Back Projection |23|. Indeed, for AT and ODT, the expression of FBP can be simplified 
into a weighting of the 1-D radial frequency measurements (e.g., as obtained with F rad ), i.e., 
a filtering of the spatial measurements, followed by a backprojection which amounts to perform 
a 2-D Fourier inversion of the weighted data extended to the whole frequency domain by zero 
padding. In AT, the frequency filtering is the so-called "ramp filter" |23|, while it can be 
shown that in ODT it corresponds to a simple Hilbert transform of the affine variations |3}|28|. 
Henceforth, due to the correspondence between LS and FBP, every solution of the ODT problem 
obtained with the LS method above will be denoted ttFBP- 



The least square method presents some strong limitations when the model (13) is severely 
ill-conditioned or when noise corrupts the observation of y. For instance, when the ODT sensing 
is associated to a weak covering of the frequency domain by the grid Vm-, i-e., for small value 
of Nq and M = NqN t < N = Nq, spurious artifacts appears in the least square estimated 
RIM ufbp- This is better understood in the simplified model where the sensing domain lives 
also on a Euclidean grid and by once again forgetting the transparent action of 0. In this case, 
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Figure 2: (left) RIM TV Model: the gradient of the RIM (represented by arrows) is non-zero only on the interface 
between the two RIM homogeneous areas no and tii. This induces a small TV norm, (right) The Shepp- Logan 
phantom, an example of the "Cartoon shape" model. 



the NDFT F reduces to a regular DFT F TCg followed by the selection of M frequencies in the N- 
dimensional Fourier domain, as represented by a restriction matrix R € {0, i} MxAr . Therefore, 
taking into account the correct normalization of the Fourier operators, we have 3> = Nq DRF rcg 
and 



y = *n = (RF I0g )*RF 



n 



4> * n, 



where is the Point Spread Function (PSF) of the frequency sampling, i.e., the convolution 
filter associated to the circulant matrix (RF rcg )* RF rcg (since R*R is diagonal). While 
is a Dirac at M = N involving ttFBP = tl, when M is significantly smaller than N, the 
PSF 4> is enlarging and contains more and more oscillations. Consequently, the solution tifbp 
deviates more and more from the pure RIM. This analysis is similar to the characterization 
of the dirty map obtained for instance in radio interferometry by aperture synthesis, which is 
another tomographic inverse problem 35 . Moreover, as the method does not take the noise 



into account, the reconstruction is generally smooth and blurred compared to the actual RIM n 



when a significant noise plagues the model (13). 



Despite these limitations, since LS/FBP is widely used in tomographic methods, we will use 
it as a standard to compare with the quality of the regularized reconstruction developed in the 
following. 



TV-^2 minimization 

In order to overcome the limitations of the LS method, we introduce a new reconstruction 
method which is both less sensitive to unwanted oscillations due to a low density frequency 



sampling Vm an d to additional observational noise rj in (11). 



In particula r, sin ce the spatial dimensions in Vm and in Cn are expected to be equal, 



i.e., 

No « N T (Sec. 3.1) we are interested in lowering the density of Vm in the Fourier plane by 
decreasing the number of angular observations Nq. In other words, with respect to this re- 
duction, we aim at developing a numerical reconstruction which makes Optical Denectometric 
Tomography 11 compressive" in a similar way other compressive imaging devices which, inspired 
by the Compressed Sensing paradigm 6,12 , reconstruct high resolution images from few (indi- 



rect) observations of its content [14,25,32,35 . This ability would lead of course to a significant 
reduction of the ODT observation time with potential impact, for instance, in fast industrial 
object quality control relying on this technology. 

This objective of compressiveness can only be reached by regularizing the ODT inverse prob- 
lem by an appropriate prior model on the configuration of the expected RIM n. Interestingly, 
the actual RIM of most human made transparent materials (e.g., optical fiber bundles, lenses, 
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optics, ..) is composed by slowly varying areas separated by sharp boundaries (material inter- 
faces) (see Fig.[2}(left)). This can be interpreted with a Bounded Variation (BV) or "Cartoon 
Shape" mod el [31| as the typical Shepp-Logan phantom in Fig. [2]-(right). Therefore, the inverse 
problem in ( |11[ ) can be regularized by promoting a small Total- Variation norm, which in its 
discrete formulation is defined as |7| 



2,1) 



where || • 1 1 2 , i is the mixed £2/^1 norm defined on any u = (1*1,1*2) £ M Arx2 as [| i*. [| 2 l = 
Ef=i(( u i)fc + ( u 2)l) 1/2 , and V : 1^ — >■ R Nx2 is the (finite difference) gradient operator. 
Reusing the 2-D coordinates of n, this operator is defined along each direction as Vn = 
(Vin, V 2 n) with (Vin) fc/ = n k+h i - n kjl and (V 2 n) fcZ = n M+ i - n k j. 

In order to obtain a reconstruction method which is also robust to additive observation noise, 



we must lighten the strict fidelity constraint implicitly used by the LS method in ( 14 ). Therefore, 
assuming the data corrupted by an additive white Gaussian noise (AWGN), we impose a data 
fidelity requirement using the ^-norm, i.e., if rj = y — <&n is known to have a bounded norm 
(or energy) \\rj\\ ^ e, we force any reconstruction candidate u to satisfy \\y — Qu\\ ^ e. The 
way e can be estimated will be explained in Sec. [6j 

Additionally to a fidelity criterion with the observed data, other requirements can be imposed 
on the reconstruction. First, we can often assume that the reference refractive index n r (i.e., 
the one of some optical fluid surrounding the object) is lower than the object RIM. Second, if 
the object is completely contained in the field-of-view Q of the ODT experiment, we can force 
any candidate RIM u to match n r on the boundary of f2, i.e., imposing u k = n r for all indices 
k belonging to the border of Cm- Up to a global value shifting the RIM (explained lates in this 
section), this is equivalent to assume that this map vanishes outside of ft. These indices are 
associated to pixels r m ^ n G Cn for which at least one of the two 2-D coordinates is equal to 
either —Nq/2 or Nq/2 — 1. The corresponding index set is denoted 9f2 for simplicity. 

Gathering all these aspects, we could propose the following reconstruction program 



tvrv-£ 2 = argmin \\u\\ T v s.t. \\y - $u\\ 2 ^ e, u y n r , ugn = n T lgn, 

ueR N 

denoting by 1 S M. N the vector of ones, i.e., the unit RIM in Cn, and recalling that vj± = Rj^v 
stands for the restriction of the components of v G ~R N to A C [N], 

However, the reconstruction can be slightly simplified by observing that the kernel of the 
sensing operator $ = &DF in ODT contains the set of constant vectors in WL N . This is a 
consequence of the vanishing affine frequency a: (which mainly defines the action of D) on 
the frequency origin, or more simply, this relies on the occurrence of the RIM gradient in the 
deflection model ([2]). 

Therefore, a change of variable u — > u — n r l does not disturb the previous reconstruction 
which can be recast as 

nxv-£ 2 = argmin ||u||tv s.t. \\y — <&w||2 ^ e, u >z 0, ugn = 0, (15) 



remembering that the true RIM estimation is actually utv-£ 2 + n rl- This last formulation has 
an interesting property. 



Lemma 5.1. If there is at least one feasible point for the constraints of (15), then the solution 
of this problem is unique. 
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Proof. Using the TV norm definition (and squaring it) , the TV-^2 optimization ( 15 ) is equivalent 
to solve 

arg min || Vu||| i s.t. \\y — «frw||2 ^ £, u y 0, uqq = 0. 

u&R N 

We remark also that the kernel of V is the set of constant vectors in M. N while the kernel of Rgn 
(defining the last constraint) is the set of vectors equal to on d£l. Therefore, Ker VnKer Rqq = 
{0}. 

Since the domain of ||V • || 1 1 is M. N , and since we assume at least one feasible point for 



the constraints of (15), we know that the program has at least one solution. Let us assume 
that there are two distinct points X\ and X2 minimizing this optimization and let us denote 
x\ = \x\ + (1 — X)x2 their convex combination for some < A < 1. By convexity, x\ is also a 
feasible point of the constraints since both minimizers X\ and X2 must satisfy the constraints. 

Moreover, Rq^x\ = RqqX2 = and x\ — X2 £ Kev Rqq. Therefore, since x\ — X2 ^ 0, 
x\ — X2 KerV and Va?i ^ Vct^- By the strict convexity of the function (/?(•) = || • W21, 
ip(Vxi) = ip(Vx2) involves 

cp{Vxi) = \(p(Vxi) + (1 - X)ip(Vx 2 ) > cp(Vx\), 

showing that x\ is a better minimizer, which is a contradiction. □ 

Therefore, we see that the uniqueness is actually reached by the stabilizing condition uqq, = 
Rdtt u = making the optimization running outside of ker V \ {0}. 

As explained in Section [7j the program (15) can be efficiently solved using proximal meth- 



ods [9] and operator splitting techniques, like the recently proposed primal-dual algorithm by 
Chambolle and Pock in |8|. 



6 NOISE IDENTIFICATION, ESTIMATION 
AND ANALYSIS 

In this section we first discuss about the different sources of noise and how to estimate the noise 
energy present in the experimental data. Then, we analyze the noise impact in both AT and 
ODT measurements. 



6.1 Noise sources 

When a real sensing scenario is being studied, such as the ODT, different sources of noise are 
present and they have to be considered when determining the global noise energy bound e in 



(15|). 

First, we have the observation noise. Under high light source intensity, the ODT mea- 
surements z produced by a Schlieren deflectometer (Fig. [lj-(bottom)) are mainly affected by 
electronic noise such as the CCD thermal noise. This induces a homoscedastic Gaussian noise 
in the measured deflection angles, i.e., with an homogeneous variance through all the measure- 
ments. By computing the 1-D Fourier transform of the ODT measurements using 0_F rac j in ([7]) 

■ 2 ) 

obs/ ' 



the corresponding noise 77 obs remains Gaussian 27 in the FDM y, i.e., r) ohs ~ M(0, a 2 



This noise, defined as the difference between the noisy FDM and the noiseless FDM y tr 
has an energy that can be bounded using the Chernoff-Hoeffding bound [18]: 



\V - ^truell 2 = ll^obsll 2 < ^obs : = tfohs ( M + C ^0' 



with high probability for c = 0(1). 
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Second, we have the modeling error that comes with every mathematical discrete repre- 
sentation of a physical continuous system. In the ODT system, this error is due to (i) the 
paraxial approximation used to formulate ([2]), (ii) the sampling of the continuous RIM and (iii) 
the discrete model itself. The modeling noise is related to the difference between the noiseless 
FDM and the sensing model <&truett = ©-D-Ftrueti.) where -Ftrue performs the exact Polar Fourier 
Transform. The modeling noise can be computed as: 

Ill/true — ^truetl|| = H^modelll < e model- 

Third, we have the NFFT interpolation noise, given by the mathematical error committed by 
estimating the polar Fourier Transform with the NFFT algorithm, i.e., the noise $ t ruei — 
To determine a bound e n st on the energy of this error we first estimate the NFFT distortion 
(i.e., without the action of D), defined as the difference between the NFFT polar Fourier 
Transform -Fapp and the true NDFT F. Theoretically, for any vector / S K . , the ^-norm of 
this distortion is bounded as \\F app f — Ff Woo ^ eC(/) = C(e||/||i), where e controls both the 
accuracy and the complexity 0(N log N/e) of the NFFT [24] (see Appendix [B]) . Assuming that 
each component of F app f — Ff is iid with a uniform distribution U([—C(f), C(/)]), and using 

the Chernoff-Hoeffding bound H20], we can estimate ||F app / - Ff\\ 2 < ^f^(M + cVM), 
with high probability for c = 0(1). Finally, e n st can be crudely computed as 

||*truen - *n|| = ||G£>(F app n - Fn)|| < |||©£)|||||F app n - Fn\\ 

_ 27r(i?r) 2 aj max eC(n) fM , „ . /77U/2 



ra vft £C(n)(M + cVM) 1/2 =: £ n fft, 

with o; max ~ \N t 5<jJ = ^ ~ representing the maximum frequency amplitude in Vm- In 
practice, because of the RIM shift n — > n — ln r explained in Sec. [5j we can bound ||n||i and 
hence C(n) with the expected RIM dynamics 5n, i.e., ||n||i ^ N5n, and we adjust e in order to 
have e n fft much lower than the other sources of noises (mainly £obs)- 

Finally, we may also have an error introduced by the instrument calibration, when deter- 
mining the exact r and 6 associated to the projections. We are going to neglect this error 
by assuming a pre-calibration process that provides an exact knowledge of these values (see 



Sec. 8.2) 



In conclusion, gathering all the previous noise identifications, we can bound the difference 
between the actual ODT measurement and the sensing model as follows: 

\\y - *n|| = \\y - y tme + y tvuc - * t ru C n + *t ruc n - *n|| 

< \\V - Vtnisll + Ill/true ~ *tru e n|| + ||*tru e n - *tl|| 

< ^obs + Emodel + Enfft = £■ 

In the following, we will neglect the error coming from the modeling £ m odel> such that only 
the measurement noise £ ODS and the NFFT interpolation noise £ n fj t are taken into account, with 
the last one considerably smaller. 

6.2 Observation noise estimation for experimental data 

When we work with experimental data, we do not have information about the variance of the 
observation noise a^ hs present in the FDM, thus we are not able to estimate the value of £obs* 
For this, we first estimate the variance a v of the noise present in the ODT measurements (z), 



denoted by t]{0,t). The Robust Median Estimator 13,34 was used for this purpose. Denoting 
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Figure 3: AT Measurement (in Fourier) (left) On the whole grid Vm for N e = 180. (right) The slice = 80° 




by ^n Q(r) = A„(t, 9) the measured signal at one specific 0, the median estimator is applied as 
follows: 

2 _ median r (r) [ 

aA ' ~ 0.6745 || h\\ 

where h{r) is a high-pass filter, * represents the discrete convolution on r and the median 
is taken for all the values of \h * ^^^(t)]. The resulting cr^(9) is the estimated variance of 
the observation noise present in the projections z. As the Fourier transform preserves the 
variance unchanged, the variance of the noise present in the FDM is the same as the one of the 
observations, i.e., <7o bs (#) = 



6.3 AT vs ODT Noise 



As we have seen in Sec. 3.3, the main difference between the AT and ODT problems is the 
appearance of the diagonal operator D in the last one. We will now analyze the impact of 
additive white Gaussian noise on the measurements regarding the presence of this operator. 

In Fig. [3] and Fig. [4] we show the magnitude of the Fourier measurements in AT and ODT, 
respectively, corresponding to the acquisition of a section of the fibers bundle (see Fig. [5]- (left)). 
For the class of images we are interested in, we can notice than in AT the magnitude presents a 
peak around lo = and then decreases significantly when the distance to the center increases, 
tending to zero in the borders (see Fig. [3]). Whereas in ODT, the presence of operator D makes 
the image intensity to be quite spread through all the pixels (see Fig. [4]). This has a direct 
impact on the reconstruction when the measurement is affected by additive Gaussian noise. As 
the noise spreads evenly through the image, the pixels that are not around uj = will be more 
affected in the AT model because their intensity is significantly lower. 
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7 NUMERICAL METHODS 



This section provides first an overview of the primal-dual algorithm we use for solving (15) 
and hence recovering the RIM from the noisy FDM. Then, the algorithm is generalized into a 
product space optimization that allows the minimization of more than two convex functions. 
Finally, we conclude the section by showing how to use this generalization to solve our ODT 
problem. 

7.1 The Primal-Dual Chambolle-Pock Algorithm 

The Chambolle-Pock (CP) algorithm (Algorithm 1 in [§]) is an efficient, robust and flexible 
algorithm that allows to solve minimization problems of the form: 

min F[Kx) + G(x), (16) 

for a linear operator K : ]R — > M. w and any variable x 6 M. N . The functions F and G belong 
to the functional sets T^IR 1 ^) and r°(IR Ar ), respectively. 

In short, CP solves the primal problem described above simultaneously with its dual problem, 
until the difference between their objective functions - the primal-dual gap - is zero. 

For any variable v G M. w , the primal-dual optimization can be formulated as the following 
saddle-point problem: 

min max (Kx, v) + G(x) — F*(v), 

where F* is the convex conjugate of function F provided by the Legendre transform F*(v) = 
ma^ m w(v,v) - F(v). 



Using the Legendre transform we obtain the primal version described in (16) and also the 
dual version as follows: 

max -F*(v) - G*(-K*v), (17) 

where K* is the exact adjoint of the operator K, such that (Kx, v) = (x, K*v). 
The CP algorithm is defined by the following iterations: 

v (k+i) _ prox^^W + vKx^), 

x (k+i) = pioXfiG ( x (k) _ M K*„(fc+i)) 5 (18) 
x (k+i) _ x (k+i) +#r x (k+i) - x (k)\ 

The quantity proxj denotes the proximal operator of a convex function / £ r°(V) for a 
certain finite dimensional vector space V [9,26]. This operator is defined as: 



proxjC := argmin/(C') + \\\C ~ C'\\ 2 , C £ V. 
C'eV 

The proximal operator admits the use of non-smooth convex functions as the TV norm, making 
the algorithm suitable to solve the TV-^2 problem described in Sec. [5j 

Most numerical methods require operator K being in a tight frame, which is not the case 
for our sensing operator The CP algorithm reduces the convergence requirements, since we 
only need to tune the step sizes \i and v such that the condition //z;|||.ff ||| 2 < 1 is true for any 
operator K. 
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7.2 Primal-Dual Gap 

The primal-dual gap, defined as the difference between the primal and dual objective func- 
tions [5]. helps us to perform a pratical analysis of the algorithm convergence. 

Setting P as the cost of the primal problem from ( 16 ) and D as the cost of the dual problem 
from (17), we can define the estimated primal-dual gap as PDG = P — D. While iterating, the 
primal objective function is greater than the dual objective function and when approching to 
the optimal solution, the difference between these objective functions tends to zero. 

If (x,v) is the output of the CP algorithm, we are interested in having the primal-dual gap 

PBG{x,v) = F(Kx) +G(x) + F*(v) + G*{-K*v), (19) 

close to zero. When the PDG is exactly zero, the strong duality holds and the solution (x,v) 
is optimal [5]. 

7.3 Product Space Optimization 

In this paper, we are interested in minimization problems containing more than two convex 
functions. In particular, we aim at solving the general optimization 

p 

min Fj(Kjx) + G(x), 

3=1 

with Kj : 1^ — > M. Wj and p + 1 the number of convex functions. Such a problem does not 
allow the direct use of the CP algorithm as described before. However, it is easy to adapt 
it by considering a p-times expanded optimization space MP N . This space is composed of 
t = (tj , ■ ■ ■ ,tp) T G W N , with tj E K . In this context, we define p — 1 bisector planes 
Hi j = {t G MP N : t\ = tj, 2 ^ j ^ p} in order to work with the following equivalent primal 
problem: 

p p 
min ^ Fj(Kjtj) + ]T (t) + ff(ti). (20) 

The CP-shape ( |16[ ) is thus recovered by working in this bigger space W N and by setting 
F{s) = YJj=i Fj(sj), with s = (sj, ■■■ , S T) T G R w= ^ w > and 8j eR w *,K = diag(Ki, • • • , K p ) G 
R W **> N and G(t) = E?= 2 *n w (*) + H(t x ). 

In this expanded optimization space, the equivalent dual problem is written (see Appendix C.l ): 

p / p 



max 

s£R w 

For the functions described above, it is easy to see that for any v > and £ = (Cf , • • • , Cp ) T £ 
M. w we have: 

fprox uF * Cx 

prox„ F „ C = 

\prox vi ,* C py 

and, for any /j, > we have (more details in Appendix D.l[ ): 

prox MG C = (i", • • • ,I JV ) T prox^ (j£,C,). 
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7.4 ODT numerical reconstruction 

Now we need to transform our TV-^ 2 problem into an expanded form in order to use the CP 
algorithm. Having two constraints, the optimization space needs to be expanded by p = 2. 



Using (20), we can reformulate the primal problem from (15) as 



min ||Vti|| 2 ,i + *c(**2) +*Po(*i) +*n 12 (t), (21) 

t=(t 1 ,t 2 )eR2JV 

where C = {v G R M : \\y - v\\ < e} and V = {u G R N : u, t ^ if i G int tt; m = if % G 317}. 
We show easily that (III]) has the shape of (|20j) with Fi(si) = ||si|| 2 ,i for si G m,Nx2 ~ ro2Af 



F 2 (s 2 ) = z c (s 2 ) for s 2 G M A/ , (t x ) = ip (*i), K x = V and K 2 = * = @DF G M MxAf . 

For building the dual problem, we need the conjugate functions of Fx, F 2 and H, which are 
easily computed using the Legendre transform. As a matter of fact, F£ is the indicator function 
onto the convex set Q = {q = (Qi,Q 2 ) G M Arx2 : ||q|| 2 ,oo ^ 1} with the mixed £oo/^ 2 -norm 



defined as ||q|| 2 ,oo = max^ J (qi)% + (q 2 )k 8 • The conjugate function F| is computed as (see 



Appendix. C.2) 



F 2 ( s 2) = *c( s 2) = (s 2 ) y + e||s 2 ||, 

while the convex conjugate of H is simply H*(s) = ^•p (— s). 
The dual optimization problem is thus defined as: 

max -iq(sx) - (s 2 ,y) - e\\s 2 \\ - zp (V*si + $*s 2 ). 

s m 2N+M 



In order to apply (18), we must compute the proximal operators of F*, F£ and H. The one 
of Fi is given by [8] 

The proximal operator of F% is determined via the proximal operator of F 2 by means of the 
conjugation property defined in [9j: 

prox^* s 2 = s 2 - ^proxi F £s 2 - 

2 v 

The proximal operator of F 2 is given by the projection onto the convex set C: 



proxi F2 s 2 = y + (s 2 - y) min (l, \\ 8 ^_ y ^ j 



The proximal operator of the function H represents a projection onto the positive orthant with 
zero borders: 

prox^£ = projp^ = 4 ' /+ . , £ G R . 

2 10 if ? G oil, 

Finally, making use of the above computations and taking $ = 1, the CP algorithm applied 



to our TV-£ 2 problem in ODT can be reduced to (see Appendix D.2) 



( S [ k+1) =prox^ («f ) + 1 /V#), 
4 fc+1) = P rox^(4 fc )+^W), 

=proj 7 , (^ fc )-^V* S f +1) + **4 /c+1) ), 1 j 
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In our experiments, the x^°\ x^ , and are initialized to zero vectors. 

In order to guarantee the convergence of the algorithm, i.e., to ensure that converges 



to the solution of (15) when k increases, we need to set \i and v such that jUi/|||.K'||| < 1. We 



are interested in giving more weight to the regularization term (TV- norm) in (15). As the 
TV-norm is contained in function F, we need to assign to the parameter v in CP a higher value 
compared to jjL. Taking v = 10 3 // provides thus more rapid convergence. The induced norm 
of the operator (|||i^|||) was computed as explained in [32] using the standard power iteration 
algorithm to calculate the largest singular value of the associated matrix K. 



The algorithm presented in ( 22 ) stops when it achieves a stable behavior, i. e., when ||aj( fe+1 ) - 



x^ ||/||a;( fc ) || ^ Th. The threshold Th is defined for ODT in the next section. In parallel, we 



analyze the convergence of the algorithm by means of the primal-dual gap as described in (19). 
If (x, s) is the output vector of the CP algorithm ( f22| ), the primal-dual gap reads 

PDG(2, a) = || V5|| 2 ,i + ic{®x) + i v (x) + iq(si) 

+ (s 2 ,y) +e||s 2 || +«p (V*si + $*s 2 ). 

Due to the presence of the indicator functions, the primal-dual gap will converge to zero 
only when their value is zero, i.e., when the argument of the indicator function belongs to the 
corresponding convex set. We will therefore evaluate a partial primal-dual gap 

pPDG(£,S) = ||V5|| 2 ,i + {s 2 ,y) + e\\s 2 \\, (23) 

and, in parallel, we monitor the value of the different indicator functions, such that we comply 
with the following conditions: 

*S G C, (24a) 
x£V , (24b) 
si G Q, (24c) 
(V*5i + $*5 2 ) G V Q . (24d) 



8 EXPERIMENTS 

In this section, the ODT reconstruction is first compared to the common tomographic (AT) 
reconstruction using the FBP method. Then, the proposed regularized reconstruction (TV-^ 2 ) 
is compared with the FBP method on synthetic and experimental ODT data. 



8.1 Synthetic Data 

Three kinds of discrete synthetic 2-D RIM are selected to test the reconstruction. They are 
defined on a 256 x 256 pixel grid (N = 256 2 ). In the first object, the RIM (n) simulates 
a 2-D section of a bundle of 10 fibers of radius 8 pixels each, immersed in an optical fluid 
(the background). The two media have a refractive index difference of 5n = 12 x 10 -3 (see 
Fig. [5]- (left)). The second object consists in a homogeneous ball centered in the pixel (154, 154) 
with a radius of 60 pixels, immersed in a liquid with 5n = 2.8 x 10~ 3 (see Fig. [5]-(right)). 
These two objects were selected in correspondence to the available experimental data we use for 
reconstruction later in this section. The third object is the well-known Shepp-Logan phantom 
(see Fig. [2J- (right)), which is a more complex image in a "Cartoon-shape" model. 

The measurements were simulated according to ( |11[ ) by means of the operator 3?, and then, 
additive white Gaussian noise ?7 obs ~ iid J\f(Q, <J 2 ob J is added in order to simulate a realistic ODT 
scenario. 
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Figure 5: Realistic refractive index map: (left) Synthetic fibers bundle and (right) Synthetic ball. 

The operator <1? is defined as *& t = &DF t , with e representing the distorsion of F e regarding 
a true operator F t me that would provide the actual NDFT. As discussed in Sec. |3.3[ the NDFT 
computational time is inversely proportional to this parameter e in O (N \og(N / e)) . Therefore, 
we need to do a compromise between an accurate and an efficient computation of the NDFT. 
For this reason we use two different operators: (i) an accurate and high dimensional operator 
<& eo = @DF eo for the acquisition, with a small eo = 10 -14 ; and (ii) a less accurate but lower 
dimensional operator 3> ei = &DF tl for the reconstruction, with ei > eo- The error caused for 
using a higher e for the reconstruction is taken into account in s D g t (see Sec. [6]). 

For each object, the ODT measurements are obtained with N T = 367 according to a varying 
number of orientations Ng, which allows to analyze the compressiveness of the reconstruction 
method. In this synthetic experiment, the orientations 6 are taken in [0, 7r) so that Ng = 360 
corresponds to two orientations per degree. Hereafter, we consider this last situation, i.e., 
59 = 7r/360 as a "full observation" scenario since, given the considered RIM resolutions, the 
discrete frequency plane is almost fully covered in this case. More generally, we say that a given 
orientation number Ng is associated to (100 of the full coverage, iVf u u being the number 

of orientations for having 56 = 7r/360. 

The reconstruction robustness with respect to the noise level has been considered by taking 
a v so that the Measurement SNR (MSNR), as measured by MSNR = 201og 10 || A||/||?7||, is 
taken in {10 dB, 20 dB, oo}. This last case with MSNR close to +oo corresponds to the noiseless 
scenario, where no Gaussian noise is added, only the NFFT interpolation error (fj^) is taken 
into account. This actually provides a high MSNR value around 270 dB. 

The reconstruction quality of n G {tifbp, titv-^I * s measured using the Reconstruction 
SNR (RSNR) measured by RSNR = 201og 10 ||n||/||n - n||. 

8.1.1 AT vs ODT 

In order to assess numerically the impact of operator D in ODT, we compare the RSNR 
between AT and ODT in similar noisy acquisition scenario. The comparison is made using 
the FBP algorithm, commonly applied in tomographic reconstructions. We analyze the impact 
of the affine frequency uj, present in ODT, via the compressiveness and noise robustness. For 
this, we focus on the reconstruction of the bundle of fibers for different number of orientations 
Ng e {4, 18, 36, 90, 180, 360}. The results are depicted in Fig. © 

In Fig. ©(left) we can see that, when no Gaussian noise is added, we obtain similar RSNR 
for both AT and ODT. The impact of the parameter uj is evident only in the convergence time, 
causing the ODT reconstruction to be 4 times slower than the AT one. However, when we add 
Gaussian noise in such a way that both data have a MSNR = 20 dB, the AT reconstruction 
presents a fast degradation while the ODT reconstruction remains almost unaffected by the 
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Figure 7: FBP vs. TV-£ 2 for different number of orientations N g with (left) MSNR = oo and (right) MSNR = 
20 dB and MSNR = 10 dB. 



noise (see Fig. [6J- (right)). These results corroborate the discussion in Sec. 6.3 



Following the discussion from Sec. 3.3, we analyze now the reconstruction of the RIM using 
a simplified ODT sensing model that is close to a classical tomographic model (12). In Fig. [6]- 
(right) we show a third curve that corresponds to the RIM reconstruction from a noisy ODT 
sensing where the diagonal operator D has been removed (ODT without D) as in (12). The 
results were obtained using the FBP algorithm and for a MSNR = 20 dB. As it was expected, 
when removing the operator D from the ODT sensing model, the reconstruction quality de- 
creases significantly compared to the results obtained with the complete ODT sensing model 
(11). Moreover, the regularized formulation TV-^2 cannot be used for this ODT reconstruction 
because the noise is then heteroscedastic. 



8.1.2 TV-^2 Reconstruction method 

The TV-^2 reconstruction is compared with the common FBP method. The reconstruction 
quality is investigated with respect to compressiveness and noise robustness. Fig. [7] presents 
comparison graphs of FBP and TV-^2 showing the RSNR vs the number of orientations Ng E 
{4, 18, 36, 90, 180, 360} for the three noise scenarios. These results correspond to the reconstruc- 
tion of the bundle of fibers for Th = 10 -5 . 

In Fig. [7]-(top left) we present the scenario without added noise, i.e., MSNR = oo. We can 
see that for a full coverage, i.e., Ng = 360, as the TV-^2 method takes into account the small 
noise coming from the NFFT interpolation error, it provides a very good reconstruction that 
outperforms by 47 dB the FBP reconstruction quality. 

The FBP method degrades rapidly when the problem is ill-posed, i.e., when the projections 
space is not fully covered, whereas the TV-^2 method maintains a high performance. By pro- 
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(a) tiFBP. 5%. RSNR = 4.78 dB. 



(b) tvrv-^. 5%. RSNR = 57.59 dB. 



x10 J x10" 




(c) n FB p. 25%. RSNR = 12.87 dB. (d) n T v-« 2 - 25%. RSNR = 64.54dB. 

Figure 8: Reconstruction images using FBP and TV-^2 reconstruction methods for MSNR = oo and different 
number of orientations. In the left column we have the FBP reconstruction results for (a) Ng = 18 and (c) 
Ng — 90. In the right column we have the TV-^2 reconstruction results for (b) Ne = 18 and (d) No — 90. 



moting a small TV- norm, the regularized method presents high compressiveness, as it can be 
observed in the graph where a high reconstruction quality is still achieved at only 5% of 360 
incident angles, obtaining a gain of 53 dB over FBP. Although the performance of the algorithm 
decreases significantly for a coverage of 1%, it still provides a higher reconstruction quality than 
FBP. 

The high compressiveness properties of the TV-^2 method are preserved when we add Gaus- 
sian noise. We are able to obtain good quality images even for a compressive and highly noisy 
sensing. With TV-l 2 , at a MSNR = 10 dB, we get a RSNR of 16 dB for a 5% radial coverage 
compared to 4dB for FBP. However, we can notice how the reconstruction quality of TV-^2 
diminishes with respect to the noiseless scenario, whereas FBP is less affected by the noise. 

Up to now, we have shown the behavior of the algorithms for different noise scenarios and 
compressiveness ratios. Let us now observe a few reconstructed images in order to appreciate 
their visual quality and the difference with the ground truth for both methods. Fig. [8] presents 
the resulting images when reconstructing the bundle of fibers in a noiseless scenario and for 
Ng = {18, 90}, which represents, respectively, a coverage of 5% and 25% of the frequency plane. 
The algorithm is set to stop when Th = 10 -5 is reached. In Table [l] we present the number of 
iterations and the time spent to obtain the results showed in Fig. [8| 

We can notice how the TV-^2 method preserves the image dynamics even for 5% of coverage, 
while FBP provides images with implausible negative values. We can also observe that some 
artifacts appear in the FBP results when the problem becomes more and more ill-posed. 

About the numerical complexities shown in Table [T] we can notice that, to reach the same 
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N e = 18(5%) 


N e = 90(25%) 


# iter 


time 


# iter 


time 


TY-£ 2 


33270 


4h31' 


10590 


lh39' 


FBP 


18930 


2h07' 


10220 


lh04' 



Table 1: Number of iterations and time for N g = {18,90} and MSNR = oo. 
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(a) n FB p. RSNR = 12.77 dB. (b) n T v-f 2 - RSNR = 26.34 dB. 
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Figure 9: Fiber reconstruction for MSNR = 20 dB and Ng = 90. Reconstructed image n using (a) FBP and (b) 
TV-^2 reconstruction methods. Difference between ground truth n and Reconstructed image n using (c) FBP 
and (d) TV-^2 reconstruction methods. 



threshold value of 10 -5 , TV-^2 needs a higher number of iterations, and therefore more time 
than FBP. However, the reconstruction quality is clearly higher when using our regularized 
method. In the case where the quality of the image reconstruction is sufficiently high, the 
threshold can be decreased to a less resctrictive value. At the end of this section we perform 
a convergence analysis for different threshold values, which allows us to choose the suitable 
threshold for a required quality or convergence time. 

In Fig. [9] we present the reconstructed images for the bundle of fibers for a moderately 
noisy sensing (MSNR = 20 dB). We also show the error images in order to provide a better 
appreciation of the difference between both methods. 

In this noisy scenario, we can notice the borders are no longer well estimated using FBP and, 
as we have coverage of only 25%, some oscillating (Gibbs) artifacts appear. On the contrary, 
the regularized method provides a good estimation on the borders, with no visible artifacts. We 
notice certain loss in the dynamics of the image, which causes a lower RSNR compared to the 
noiseless scenario. 

Let us show now some results obtained for the other two synthetic images. Fig. 10 and 
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(a) n FB p. RSNR = 22.03 dB. 



(b) xiTV-ti- RSNR = 28.83 dB. 



xKT* x10" 




(c) |n-np B p|. (d) |n-tvrv-fel- 

Figure 10: Sphere reconstruction for MSNR = 20 dB and Ne = 90. Reconstructed image n using (a) FBP and 
(b) TV-^2 reconstruction methods. Difference between ground truth n and Reconstructed image n using (c) FBP 
and (d) TV-I2 reconstruction methods. 



Fig. 11 present the results obtained using FBP and on the sphere and the Shepp-Logan 

phantoms, respectively. This comparison was performed for MSNR = 20 dB and for Ng = 90, 
i.e., a coverage of 25% of the frequency plane. 

The regularized method provides a better image dynamics for these phantoms than when 
reconstructing the fibers for the same noise level. For these phantoms, it can also be observed 
that FBP reconstructions present a poor estimation on the borders and oscillating artifacts. 
Moreover, the error image shows a higher discordance with respect to the actual image. 

Table [2] presents a more complete comparison of the RSNR obtained using FBP and TV-^2 
on the three synthetic images. The methods are analyzed for three scenarios, one noiseless with 
MSNR = 00, and the other two with noise such that we have MSNR = 20 dB and MSNR = 
10 dB. Results are presented for Th = 10~ 5 and N$ = 90, i.e., a coverage of 25% of the frequency 
plane. 





RSNR[dB] 


MSNR = 00 


MSNR = 20 dB 


MSNR = 10 dB 


TV-h 


FBP 


TV-£ 2 


FBP 


TV-£ 2 


FBP 


Fibers 


64.54 


12.87 


26.34 


12.77 


19.20 


11.91 


Sphere 


34.65 


22.34 


28.83 


22.03 


21.85 


19.61 


Shepp-Logan 


39.13 


14.19 


35.38 


14.04 


25.33 


12.87 



Table 2: Comparison of the different RSNR obtained using FBP and TV- 
N = 90. 



on the three synthetic images for 
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(a) n FB p. RSNR = 14.04 dB. 



(b) n 



TV-l 2 



RSNR = 35.38 dB. 
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Figure 11: Shepp-Logan phantom reconstruction for MSNR = 20 dB and Ne = 90. Reconstructed image n using 
(a) FBP and (b) TV-^2 reconstruction methods. Difference between ground truth n and Reconstructed image n 
using (c) FBP and (d) TV-^2 reconstruction methods. 
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Figure 12: Convergence results when reconstructing the bundle of fibers with Ng = 360 and MSNR = 20 dB. The 
progress of (top left) the threshold value Th and (top right) the RSNR, along the iterations, (bottom left) The 
PDG (i) and (bottom right) the condition (ii) of the primal-dual gap corresponding to the indicator function ic ■ 



When comparing the behavior of the algorithm for the different synthetic images, we can 
notice that TV-^2 method outperforms the FBP method for all cases. When using the TV-^2 
method it is important to see that, for the noiseless scenario, the fibers are best reconstructed, 
preserving the image dynamics; while for the noisy scenarios, the Shepp-Logan led to a better 
image reconstruction. 

8.1.3 Algorithm Convergence 

Finally, we analyze the convergence of the algorithm by studying the evolution of the primal- 
dual gap and of the RSNR for different threshold values. For this, we use the reconstruction of 
the bundle of fibers for N e = 360 and MSNR = 20 dB. 



The evolution of ||a:( fe+1 ) — \\/\\x^ || along the iterations is depicted in Fig. 12 -(top left) 
and it helps us to analyze the stability of the algorithm. We can see the curve is not smooth, 
which indicates a non stable behavior mainly due to a bad conditioning of the global operator 
K in the product space optimization. This could be improved by a preconditioning procedure 
as described in [2|[29] . Fig. [l2|-(top right) presents the RSNR evolution along the iterations. 
The evolution of the primal-dual gap along the iterations on the current primal and dual 



solution is also analyzed by evaluating the partial PDG (pPDG) from (23) and by monitoring 
the compliance of the conditions required by the indicator functions (|24[). The behavior of the 



partial PDG can be observed in Fig. 12 (bottom left). Fig. 12 -(bottom right) presents the value 



of || y — «&sb|| compared to the value of e, which allows us to analyze the compliance of the 



indicator function %c in condition (24a). The other conditions are not presented here but we 
observed that they are satisfied in few iterations. 

Table [3] presents the values of some convergence parameters (Time, RSNR, pPDG, ...) for 
different threshold values. The parameter pPDG % refers to the relative pPDG, computed as 
pPDG% = 100((P — D)/(P + D)). These results show that a lower threshold provides higher 
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reconstruction quality but significantly increases the number of CP iterations. In the meantime, 
the data consistency \\y — «&sc|| — e also decreases. In this specific reconstruction, setting the 
threshold to O(10" 5 ) or running more than 1000 CP iterations guarantees a RSNR higher than 
31 dB. 



Th 


# iter 


Time 


RSNR [dB] 


pPDG 


pPDG [%] 


\\y — $>x\\ — e 


10- 4 


400 


8' 


29.52 


0.64 


4.43 


3.7 x 10" 3 


10~ 5 


900 


18' 


31.07 


0.31 


2.19 


2.2 x 10~ 3 


10" b 


5630 


2h 


35.4 


0.04 


0.29 


6.9 x 10~ 4 


10- 7 


32720 


llh 


40.26 


0.005 


0.04 


2.1 x 10~ 4 



Table 3: Convergence results when reconstructing the 10 fibers synthetic image with Ng — 360 and MSNR = 
20 dB. 



8.1.4 TV-^2 problem with one constraint 



We analyze now the impact of the constraints added to the TV-^2 formulation in (15). For 
this, we solve a TV-^2 problem without the positivity and zero borders constraints, which were 
proved to guarantee the existence of an unique solution. We call this problem the incomplete 
program and we compare it with the complete TV-^2 program from Sec. [5j For the 
noiseless scenario, we obtain better reconstruction quality when using the complete version. 
For the noisy scenarios, the reconstruction results are similar for both versions of the TV-^2 
program. However, the reconstruction for the incomplete program is slower, taking 5 times more 
iterations to converge than what the complete program takes. We conclude that, in addition to 
guaranteeing the unicity of the solution, the extra constraints stabilize also the convergence of 
the reconstruction. 



8.2 Experimental Data 

The reconstruction algorithm was tested with two particular transparent objects similar to the 
synthetic data studied in the previous section: a bundle of 10 fibers and a homogeneous sphere, 
both immersed in an optical fluid. The reconstruction is based on N T = 696 parallel and equally 
spaced light rays. The experimental setup is based on the Schlieren Deflectometric Tomography 
described in Sec. [21 

A 696 x 523 pixels CCD camera was used for the acquisition, covering a field of view of 
3.25mm x 2.43mm. This corresponds to N T = 696 parallel light rays and 523 2-D slices, which 
leads to 5t = 4.7 x 10 _3 mm, and thus to 5r = 4.7 x 10 _3 mm. 

The experimental configuration leads to a calibration problem. As the object is rotating, 
the center of rotation may change within a small range. Therefore, a calibration method was 
implemented as a post-acquisition phase to determine where the projected image is actually 
centered. The method computes the location of the maximum and minimum deflection values 
in the direction of r for all projection angles. Then, the center of the projection image is given 
by the mean value between these locations. 

The characteristics of the objects and the reconstruction results obtained from the measured 
data are given below. 



Bundle of fibers 

The first measured object is a bundle of 10 fibers of 200 fim diameter each. The refractive 
index difference with respect to the solution where the fibers are immersed is 5n = 12 x 10 -3 . 
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The experimental data was measured for 60 angular positions over 360 degrees (i.e., 17% of 
measurements). Fig. 13 shows the reconstruction results obtained using both FBP and 
reconstruction methods. 



x10" 3 x10" 




x(mm) x(mm) 

Figure 13: Reconstructed bundle of 10 fibers for 60 angular positions. 2D distribution using (top left) FBP and 
(top right) TV-4. ID profile along y = 1.218 mm using (bottom left) FBP and (bottom right) TV-£ 2 . 

We observe a similar behavior to the one found in the synthetic reconstruction. We have a 
more accurate image dynamics recovery using our regularized method, whereas with FBP the 
reconstruction presents several artifacts and we observe the presence of implausible negative 
values. It is important to notice that the preservation of the image dynamics depends mainly 
on the noise estimation and on the proper definition of the constants included in the operator D 
(see Sec. [3]). When considering the appropriate constants, we are able to make an equivalence 
between the physical problem and its discrete mathematical formulation. 

Note that the interior of the reconstructed object is not entirely homogeneous, such recon- 
struction errors are still present because the modeling error was not taken into consideration. 
Given a high refractive index difference 5n, the first order paraxial approximation is no longer 
valid and thus the modeling error becomes more important in 0. The actual light ray tra- 
jectory could be estimated and inserted in an iterative process as done by Antoine et al. HI. 
However, we would then loose the fast computation of the forward model in the frequency do- 
main. Moreover, a section of the bundle of fibers represents a complex image to reconstruct, as 
the light enters and comes out of multiple fibers causing an error propagation. 

As the dimensions of the problem are higher in the experimental setup, the reconstruction 
time increases. For a threshold of 10~ 6 , FBP converges in 50 x 10 3 iterations and TV-^2 in 
44 x 10 3 iterations. 

Homogeneous Sphere 

The second observed object consists of a homogeneous sphere with a diameter of 1.55 mm. The 
difference of refractive index between the sphere and the optical fluid where it is immersed is 
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Sn = 2.8 x 1(T 3 . The deviation map was measured for Nq = 45 angular positions over 360 
degrees (i.e., 13% of measurements). Fig. 14 shows the reconstruction results obtained when 
using both FBP and TV-^2 reconstruction methods. 




We note similar results as the obtained when reconstructing the fibers: a more accurate 
image dynamics recovery using our regularized method than FBP. However, when reconstructing 
the sphere with the TV-^2 method, we can notice that the interior of the recovered object is 
more homogeneous and closer to the expected value. This is due to the lower refractive index 
difference between the sphere and the liquid where it is immerged, which provides a smaller 
deviation angle and, therefore, a smaller modeling error. Moreover, the sphere represents a less 
complex image compared to the fibers, as the light enters and comes out of only one object, 
thus there is no error propagation. 

The results are shown for a threshold Th = 10 -6 , where FBP converges in 50 x 10 3 iterations 
and TV-^2 in 40 x 10 3 iterations. 



9 CONCLUSIONS 

We have demonstrated how regularized reconstruction methods, such as TV-^2> can be used in 
the framework of Optical Deflectometric Tomography in order to tackle the lack of deflectometric 
observations and to make the ODT imaging process robust to noise. The proposed constrained 
optimization problem shows significant improvements in the reconstructions, compared to the 
well-known Filtered Back Projection. The results confirm that when dealing with a compressive 
setting, the Total Variation regularization and the prior information constraints (positivity and 
FoV restriction) help in providing a unique and accurate estimation of the RIM, promoting sharp 
edges and preserving the image dynamics. By working with the Chambolle-Pock algorithm we 
exploit the advantages of proximal operators and of primal-dual algorithms, and their flexibility 
to integrate multiple constraints. We have also shown that the use of the fast NFFT algorithm 



30 



efficiently approximates (with a controlled error) the ODT sensing model involving the polar 



Noticeably, there still exist some artifacts in the experimental data reconstruction, coming 
from the modeling error. In order to handle this problem, we should no longer assume a linear 
light propagation. The actual curved light trajectories, depending themselves on the RIM 
through the eikonal equation, must be traced. However, such a model improvement breaks the 
fast computability of the forward sensing operator through the Fourier domain and knowing 
how to insert this new scheme in our reconstruction method is a matter of future works. 

We have also noticed the long convergence time of the algorithm. As the CP algorithm has 
been proved to work slowly when the problem is badly conditioned 1 2 , 29 1 , convergence results 
could be improved by preconditioning our global operator K. We note also that most of the 
algorithms are implemented in Matlab, except for the NFFT algorithm which is compiled in 
C++. This makes the implementation slower, and could be improved for instance by the use 
of parallel processing techniques. 

We finally remark that the Optical Deflectometric framework treated in this paper can be 
also applied to other imaging techniques, such as the X-ray phase contrast tomography |28| . 
With the use of the energetic X-ray light, the first order paraxial approximation involving linear 
light trajectory is no longer needed and the proposed methods are expected to provide better 
results. 
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A DEFLECTOMETRIC FOURIER SLICE THEOREM 



NDFT. 



Using the notations of Sec. 2.1, we must prove that 



y(u>,9) 



n r 



n(cjp e ). 



By definition of y in ([3]), we have 




Vn(r) 




— 2iri ra j2 



d rdr 




However, for any function / : K — > C integrable on M 




Therefore, we compute easily for any a, £ G M 2 




Setting a = p e and £ = ojp e , we find finally 
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B NON-EQUISPACED FOURIER TRANSFORM 
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The non-equispaced Fourier Transform (NFFT) allows us to compute rapidly, i.e., in 0(N log N) , 
the NDFT defined in Q of a function defined on Cn. This computation is performed with a 
controllable error which can be further reduced by increasing the computational time. 

In a nutshell, the NFFT algorithm replaces the equivalent matrix multiplication / = Ff 
of |9| by s = Ff, where F has fast matrix-vector multiplication computation. In this scheme, 
the discrepancy between / and s, as measured by E oc (f) := \\f — s||oq, is controlled and kept 
small. 

More precisely, the matrix F starts by embedding ~R N in a bigger regular space M. n of size 
n = ng, with no > Nq. This is obtained from the multiplication of / with a matrix W G M nxAr 
that performs a weighting of / by a vector w G (component wise) followed by a symmetric 
zero-padding on each side of the function domain Cn- Once in M. n , the common DFT matrix 
F n is applied. It can be computed with the FFT algorithm in order to obtained an oversampled 
Fourier transform of / G R . Finally, a sparse matrix V G M. Mxn multiplies the output of F n 
in order to end in the M dimensional space V . Each row of V corresponds to the translation 
in the 2-D Fourier domain of a compact and separable 2-D filter ip on one specific point of the 
non-regular grid Pm- As a final result, the matrix F is thus factorized in 

F = VF n W. 

Without entering into unnecessary technicalities, the NFFT scheme is characterized by a 
precise connection between the function ip defining V and the weighting performed by W. In 
particular, each component of w is actually set to the inverse of the Fourier transform of a filter 
ip, while ip is a periodization of (a truncation of) the same filter. 

There exist several choices of windows ip/ip associated to different numerical properties {e.g., 
localized support in frequency and time, simple precomputations of the windows, ...). We select 
here the translates of Gaussian bells [33], involving a Gaussian behavior for (p{k), which provides 
fast error decay for E oc (f). 

In particular, denoting by a = n/N ^ 1 the oversampling factor and using the FFT for 
matrix- vector multiplications involving F n , the total complexity T of the multiplications of F 
or F with vectors is 

T{F) = 0(nlog n + N(^Ek) log 1/e), 

if we impose E OQ (f) ^ 4||/||i e. For a fixed a > 1, this reduces to T(F) = O(NlogNfe), which 
is far less than a direct computation of the DFT in O(MN) computation even for small value 
of e. 



C CONVEX CONJUGATE FUNCTIONS 

C.l Convex conjugate of the function G. 

Given a vector t = (t^,^)^ G M? N , we can define the function G{t) = H(t) + 2(^1^2), with 
the bisector plane 111,2 = {t : t± = t2}. Therefore, for a vector u = (u 1 [,u r ^) T G M. 2N , the dual 
function G*(t) can be computed as follows: 

G*(t) = max (u,t) - G(t) 

= max (iii, t\ + £2) — H(ux) 

u: 1*1=1x2 

= (;~);r( tl + t2 ). 
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C. 2 Convex Conjugate of the indicator function % c . 

Given the indicator function ic{u) of the convex set C = {v G R M : \\v — y\\ ^ e}, its dual 
function can be computed via the Legendre transform as follows: 

i*c{v) = max (v, u) — ic(u) 

= {v,y)+ max {v,u — y) 

u:\\u—y\\^e 

= (v,y) + max (v,b). 

b:\\b\\^e 

The value of {b : \\b\\ ^ e} that maximizes the last expression is b = and we get i* c {v) = 

(v,y) +e\\v\\. 

D DETAILS ON THE PRODUCT SPACE OPTIMIZATION 

D. l Proximal operator of the function G 

For every t = (ti, ■ ■ ■ , t p ) G R pN , the function G(t) is defined as: 

G(t) = J> ni ,,(t)+#(ti), 

i=2 

where, for j£ [p], Ili j denotes the bisector plane Hij = {t G M pAr : t\ = tj}. 
Then, for £ = (£ l5 • • • , £ p ) G M pAr , the proximal operator of G reduces to 

t = prox MG £ = argmin fiH(t\) + ^HC - *|| 2 - 

t: ti = ---=tp 



u = argmin p,H(u) + | ||£ • — u\ 



As all the tj are equal, we necessarily have t = (I N , • • • ,I Ar ) T u with it G 1^ defined by 

,IH,j-W 2 

= argmin + § [||Cil| 2 + ■ ■ ■ + ||C P I| 2 - 2pl T u + p\\u\\ 2 ] 

ueR N 

= arg min fj,H (u) + | [p||CI| 2 — 2pC, T u + p||ii|| 2 ] 

= P rox ( M / P )// C- 

with £ = I £V ^ G 1^ and where we used the fact that we can always subtract or add constants 
in the minimization without disturbing its solution. Denoting by 1^ = (1^, • • • ,1^) G M. NxpN 
the operator such that lp £ = pC, this provides finally the compact notation 

N\T 1,7V, 



prox MG C = (ip Y prox^ 
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D.2 Formulation of Chambolle-Pock algorithm in OD 



We take the CP algorithm as described in Eq. ( 18 ) and we translate it into our OD problem in 
the product space, to obtain the following: 



ffc+1) 



(fc+1)^ 



prox^i(tS fc) -/,V* S f +1) +4 fc) 



M*4 fe+1) ), 



We have the function i2"(t) = 2-p (t) an d its proximal operator proxn# £ = proj-p () C- As ti = t 
we can relabel the variable as = t[ k ^ 



2: 



4 fc) and a;( fc+1 ) 



ti = <2, thus we can relabel the variable as x^ 
obtain the following algorithm: 



i(fc) 



H — c 5 



In the same way, 



t ( 2 k) and 



r(fe+l) _ ?(*+!) 



We 



,(*+!) 

^(fc+i) 



prox^»(4 A) 



proj Po (ajW - 
2a; (fe+i) _ x (k) 



(V*sf +1) + **4 fc+1) )), 
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